In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.
library(dplyr) # Data Manipulation
Warning: package ‘dplyr’ was built under R version 4.1.1
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal,
union
library(ggplot2) # 2-D Plot
Warning: package ‘ggplot2’ was built under R version 4.1.1
library(plotly) # 3-D Plot
Warning: package ‘plotly’ was built under R version 4.1.1
Attaching package: ‘plotly’
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
library(ggfortify) # used for biplots of PCAs
Warning: package ‘ggfortify’ was built under R version 4.1.1
library(vegan) # Ecological analysis mega-package
Warning: package ‘vegan’ was built under R version 4.1.1
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.1
Loading required package: lattice
This is vegan 2.5-7
library(RMTRCode2)
# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
obj <- deparse(substitute(obj))
if(exists(obj, envir = env)) {
rm(list = obj, envir = env)
}
}
First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)
In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 5 did not.
The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.
Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.
To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Every vertical line is a species introduction or extinction by neutral dynamics.
Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Mean = mean(Richness),
SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
Species, collapse = ", "), split = ", ", fixed = TRUE))))),
Gamma = unlist(lapply(strsplit(
SpeciesTotal, split = ", ", fixed = TRUE), length))
) %>% tidyr::pivot_longer(
cols = c(Mean, Gamma),
names_to = "Aggregation",
values_to = "Richness"
),
mapping = ggplot2::aes(
x = Time,
y = Richness,
linetype = Aggregation
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value:
mean(Diversity$Richness)
[1] 4.458893
(If we consider the first 10,000 time units as burn-in, we instead see a value of
mean(Diversity$Richness[Diversity$Time > 10000])
[1] 4.587217
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 8043 row(s) containing missing values (geom_path).
Warning: Removed 2297 row(s) containing missing values (geom_path).
Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.
# ggplot2::ggplot(
# Diversity %>% dplyr::filter(Environment < 4),
# ggplot2::aes(
# x = Richness,
# y = Entropy,
# group = Environment,
# color = Time
# )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
# alpha = "none"
# )
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 10986 row(s) containing missing values (geom_path).
Warning: Removed 2831 row(s) containing missing values (geom_path).
Evenness helps highlight that this is a high variance process but with a relatively constrained mean.
We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only
unique(EnvDiversity$Species[EnvDiversity$Richness > 5])
[1] 12 14 22 30
are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7
1 65970 14384 439 0 0 0 0 0
2 34776 31378 11352 3287 0 0 0 0
3 63206 13000 1767 2820 0 0 0 0
4 46367 21078 11443 1617 288 0 0 0
5 30648 31553 16542 1341 709 0 0 0
6 72061 8732 0 0 0 0 0 0
7 75522 5271 0 0 0 0 0 0
8 64759 15508 526 0 0 0 0 0
9 77919 2874 0 0 0 0 0 0
10 17817 38679 20192 4105 0 0 0 0
11 52666 24892 3235 0 0 0 0 0
12 9669 2676 9781 26641 18597 9922 2294 1213
13 64640 15703 450 0 0 0 0 0
14 12094 21314 14556 16445 15148 931 305 0
15 68564 10732 1497 0 0 0 0 0
16 69690 10903 200 0 0 0 0 0
17 30636 33714 15505 938 0 0 0 0
18 69214 11322 257 0 0 0 0 0
19 62460 17599 734 0 0 0 0 0
20 74867 3685 1849 392 0 0 0 0
21 47906 25702 4185 3000 0 0 0 0
22 3180 13477 20151 21766 16740 5305 174 0
23 65137 14940 716 0 0 0 0 0
24 65792 12313 2688 0 0 0 0 0
25 66085 11109 3599 0 0 0 0 0
26 41654 19402 19201 536 0 0 0 0
27 48808 20447 10279 1259 0 0 0 0
28 68983 11810 0 0 0 0 0 0
29 34631 26993 17044 2125 0 0 0 0
30 5561 19065 20682 11364 11597 7963 3301 1260
31 65019 15665 109 0 0 0 0 0
32 19327 26015 29600 5851 0 0 0 0
33 15128 23405 24345 10059 7856 0 0 0
34 76432 4361 0 0 0 0 0 0
35 73829 6080 884 0 0 0 0 0
36 71050 9743 0 0 0 0 0 0
37 38199 22369 15663 4486 76 0 0 0
38 28752 35327 16714 0 0 0 0 0
39 77008 3785 0 0 0 0 0 0
40 49953 27129 3711 0 0 0 0 0
41 79838 955 0 0 0 0 0 0
42 78312 2481 0 0 0 0 0 0
43 32813 23892 17056 6661 371 0 0 0
44 39700 23980 15554 1559 0 0 0 0
45 80793 0 0 0 0 0 0 0
46 48461 27669 4663 0 0 0 0 0
47 12053 36783 14826 15182 1949 0 0 0
48 63151 16980 662 0 0 0 0 0
49 73125 7668 0 0 0 0 0 0
50 79860 933 0 0 0 0 0 0
51 80255 538 0 0 0 0 0 0
52 73946 6847 0 0 0 0 0 0
53 77581 3212 0 0 0 0 0 0
54 76399 4394 0 0 0 0 0 0
55 55217 23951 1625 0 0 0 0 0
56 72564 8229 0 0 0 0 0 0
57 57443 22644 706 0 0 0 0 0
58 78988 1805 0 0 0 0 0 0
59 35315 37942 7536 0 0 0 0 0
60 67881 12912 0 0 0 0 0 0
61 52226 24749 3818 0 0 0 0 0
62 62517 16663 1613 0 0 0 0 0
63 46937 31441 2415 0 0 0 0 0
64 68974 11213 606 0 0 0 0 0
65 80793 0 0 0 0 0 0 0
66 67130 12561 1102 0 0 0 0 0
67 37887 27036 11834 4036 0 0 0 0
68 51139 27047 2205 402 0 0 0 0
69 61791 16492 2510 0 0 0 0 0
70 73231 7354 208 0 0 0 0 0
71 22117 22512 20026 14172 1966 0 0 0
72 10060 26081 33450 7590 3612 0 0 0
73 73696 7097 0 0 0 0 0 0
74 80793 0 0 0 0 0 0 0
75 77528 3008 257 0 0 0 0 0
76 34905 34766 10606 516 0 0 0 0
77 54799 21256 4116 622 0 0 0 0
78 57557 21260 1976 0 0 0 0 0
79 46966 28430 5397 0 0 0 0 0
80 46394 30478 3921 0 0 0 0 0
81 80793 0 0 0 0 0 0 0
82 72153 8640 0 0 0 0 0 0
83 57732 21901 1160 0 0 0 0 0
84 73061 7732 0 0 0 0 0 0
85 78108 2685 0 0 0 0 0 0
86 64841 15952 0 0 0 0 0 0
87 75894 4488 411 0 0 0 0 0
88 77088 3705 0 0 0 0 0 0
89 48328 32465 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0
91 36647 30720 11100 1868 458 0 0 0
92 54513 24626 1654 0 0 0 0 0
93 73946 6847 0 0 0 0 0 0
94 69274 10406 1113 0 0 0 0 0
95 43213 20293 14643 2644 0 0 0 0
96 65015 15778 0 0 0 0 0 0
97 75428 5365 0 0 0 0 0 0
98 65968 13004 1821 0 0 0 0 0
99 78629 2164 0 0 0 0 0 0
100 52405 15982 9604 2802 0 0 0 0
(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.03749 0.03273 0.03095 0.02751 0.02501 0.02482
The amount of explained variation is as follows.
sum(summary(PCA)$importance[2, ])
[1] 0.99994
The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
We can try again with presence-absence data instead.
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.05537 0.04584 0.04483 0.03901 0.03578 0.03410
So it is a bit better, but not by much.
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover
sum(summary(PCA)$importance[2, ])
[1] 0.99994
of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.
To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)
My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.
We can try again with presence absence instead.
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)
ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)
We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.
load(file.path(
"..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Every vertical line is a species introduction or extinction by neutral dynamics. Note that, due to dispersal, these should lose some of their abundance to spreading, but if that spreading is faster than the population can recoup its losses, than it can invade and be wiped out (since invadability does not consider spatial dynamics, only local ecological dynamics). The primary notes then are that there are a large number more events appearing to occur (due to the spatial arrangement), but the system as a whole is stabler (since the species that get knocked out can be recaptured from space as well).
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Mean = mean(Richness),
SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
Species, collapse = ", "), split = ", ", fixed = TRUE))))),
Gamma = unlist(lapply(strsplit(
SpeciesTotal, split = ", ", fixed = TRUE), length))
) %>% tidyr::pivot_longer(
cols = c(Mean, Gamma),
names_to = "Aggregation",
values_to = "Richness"
),
mapping = ggplot2::aes(
x = Time,
y = Richness,
linetype = Aggregation
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1123 row(s) containing missing values (geom_path).
Warning: Removed 110 row(s) containing missing values (geom_path).
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7 8 9 10
1 80793 0 0 0 0 0 0 0 0 0 0
2 80624 0 1 1 2 1 1 2 0 2 159
3 80793 0 0 0 0 0 0 0 0 0 0
4 80047 0 0 2 0 1 0 0 1 0 742
5 79905 0 6 0 6 2 5 1 0 11 857
6 80793 0 0 0 0 0 0 0 0 0 0
7 80793 0 0 0 0 0 0 0 0 0 0
8 80793 0 0 0 0 0 0 0 0 0 0
9 80793 0 0 0 0 0 0 0 0 0 0
10 80793 0 0 0 0 0 0 0 0 0 0
11 80793 0 0 0 0 0 0 0 0 0 0
12 680 0 0 0 0 0 0 0 0 1 80112
13 80639 1 1 1 1 1 1 2 1 1 144
14 84 0 0 0 0 0 0 0 0 0 80709
15 80793 0 0 0 0 0 0 0 0 0 0
16 80793 0 0 0 0 0 0 0 0 0 0
17 80793 0 0 0 0 0 0 0 0 0 0
18 80793 0 0 0 0 0 0 0 0 0 0
19 80793 0 0 0 0 0 0 0 0 0 0
20 80793 0 0 0 0 0 0 0 0 0 0
21 80793 0 0 0 0 0 0 0 0 0 0
22 62 0 0 0 0 0 0 0 0 1 80730
23 80793 0 0 0 0 0 0 0 0 0 0
24 80546 0 1 0 0 0 0 1 0 0 245
25 80793 0 0 0 0 0 0 0 0 0 0
26 80793 0 0 0 0 0 0 0 0 0 0
27 80793 0 0 0 0 0 0 0 0 0 0
28 80793 0 0 0 0 0 0 0 0 0 0
29 80721 1 0 0 1 0 0 0 0 1 69
30 69391 9 2 8 1 8 5 2 8 12 11347
31 80793 0 0 0 0 0 0 0 0 0 0
32 500 0 0 0 0 0 1 0 0 1 80291
33 66491 5 15 10 14 2 14 5 14 17 14206
34 80793 0 0 0 0 0 0 0 0 0 0
35 80332 1 0 1 1 1 1 1 0 2 453
36 80793 0 0 0 0 0 0 0 0 0 0
37 46279 1 1 2 1 1 3 1 2 4 34498
38 80793 0 0 0 0 0 0 0 0 0 0
39 80793 0 0 0 0 0 0 0 0 0 0
40 80793 0 0 0 0 0 0 0 0 0 0
41 80793 0 0 0 0 0 0 0 0 0 0
42 80793 0 0 0 0 0 0 0 0 0 0
43 27630 2 2 0 4 0 0 1 4 7 53143
44 60232 1 8 3 3 7 2 7 5 14 20511
45 80793 0 0 0 0 0 0 0 0 0 0
46 80793 0 0 0 0 0 0 0 0 0 0
47 52191 5 8 3 7 3 7 5 9 12 28543
48 80793 0 0 0 0 0 0 0 0 0 0
49 80793 0 0 0 0 0 0 0 0 0 0
50 80793 0 0 0 0 0 0 0 0 0 0
51 80793 0 0 0 0 0 0 0 0 0 0
52 80793 0 0 0 0 0 0 0 0 0 0
53 80793 0 0 0 0 0 0 0 0 0 0
54 80793 0 0 0 0 0 0 0 0 0 0
55 78735 1 1 2 0 1 2 0 0 2 2049
56 80793 0 0 0 0 0 0 0 0 0 0
57 79020 0 0 1 0 0 0 0 1 0 1771
58 80793 0 0 0 0 0 0 0 0 0 0
59 80793 0 0 0 0 0 0 0 0 0 0
60 80793 0 0 0 0 0 0 0 0 0 0
61 80466 0 1 0 0 1 1 1 0 1 322
62 80793 0 0 0 0 0 0 0 0 0 0
63 80793 0 0 0 0 0 0 0 0 0 0
64 80793 0 0 0 0 0 0 0 0 0 0
65 80793 0 0 0 0 0 0 0 0 0 0
66 80793 0 0 0 0 0 0 0 0 0 0
67 80793 0 0 0 0 0 0 0 0 0 0
68 79023 1 5 4 7 5 6 4 4 4 1730
69 78135 2 2 1 1 1 6 1 2 5 2637
70 80793 0 0 0 0 0 0 0 0 0 0
71 1394 0 0 0 0 0 0 0 0 1 79398
72 292 1 0 0 0 0 1 0 1 0 80498
73 80793 0 0 0 0 0 0 0 0 0 0
74 80793 0 0 0 0 0 0 0 0 0 0
75 80793 0 0 0 0 0 0 0 0 0 0
76 78676 1 4 1 4 0 6 1 3 6 2091
77 79765 0 0 1 0 0 1 0 0 1 1025
78 79253 0 0 1 0 1 0 0 0 0 1538
79 80793 0 0 0 0 0 0 0 0 0 0
80 80793 0 0 0 0 0 0 0 0 0 0
81 80793 0 0 0 0 0 0 0 0 0 0
82 80793 0 0 0 0 0 0 0 0 0 0
83 80793 0 0 0 0 0 0 0 0 0 0
84 80793 0 0 0 0 0 0 0 0 0 0
85 80793 0 0 0 0 0 0 0 0 0 0
86 64862 3 6 4 7 4 6 3 8 6 15884
87 80793 0 0 0 0 0 0 0 0 0 0
88 80793 0 0 0 0 0 0 0 0 0 0
89 80793 0 0 0 0 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0 0 0 0
[ reached getOption("max.print") -- omitted 10 rows ]
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.21750 0.12185 0.08092 0.05943 0.05215 0.04468
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
sum(summary(PCA)$importance[2, ])
[1] 0.99996
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)
ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.
load(file.path(
"..", "experiments", "MNA-FirstAttempt1e+06-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[1:(length(toRemove)/2)][toRemove[1:(length(toRemove)/2)]] <- 0
result$Abundance[(length(toRemove)/2 + 1):length(toRemove)][toRemove[(length(toRemove)/2 + 1):length(toRemove)]] <- 0
ifrm(toRemove)
With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.
LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1,
to = nrow(result$Abundance),
by = 10), c(1, 2:101)]) -> obj;
obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(ifelse(x != 0, x * log(x), 0))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Species = species,
Environment = i)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
Evenness = Entropy / log(Richness)
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Richness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Mean = mean(Richness),
SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
Species, collapse = ", "), split = ", ", fixed = TRUE))))),
Gamma = unlist(lapply(strsplit(
SpeciesTotal, split = ", ", fixed = TRUE), length))
) %>% tidyr::pivot_longer(
cols = c(Mean, Gamma),
names_to = "Aggregation",
values_to = "Richness"
),
mapping = ggplot2::aes(
x = Time,
y = Richness,
linetype = Aggregation
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Entropy,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Entropy = mean(Entropy)
),
mapping = ggplot2::aes(
x = Time,
y = Entropy
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 1855 row(s) containing missing values (geom_path).
Warning: Removed 259 row(s) containing missing values (geom_path).
plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
mode = "lines", opacity = 1, line = list(color = ~Time))
ggplot2::ggplot(
Diversity,
ggplot2::aes(
x = Time,
y = Evenness,
color = factor(Environment),
alpha = ifelse(Environment == 1, 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Evenness = mean(Evenness)
),
mapping = ggplot2::aes(
x = Time,
y = Evenness
),
color = "black",
inherit.aes = FALSE
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
)
Warning: Removed 2594 row(s) containing missing values (geom_path).
Warning: Removed 329 row(s) containing missing values (geom_path).
EnvDiversity <- lapply(
1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
environments <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
data.frame(Time = time,
Richness = richness,
Abundance = abundSum,
Species = i,
Environments = environments)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
EnvDiversity <- dplyr::bind_rows(EnvDiversity)
ggplot2::ggplot(
EnvDiversity %>% dplyr::filter(Richness > 1),
aes(x = Time, y = Richness, color = Species)
) + geom_point(
alpha = 0.01, size = 3
) + guides(
color = "none"
)
with(EnvDiversity %>% mutate(
Time = floor(Time)
) %>% group_by(
Time, Species
) %>% summarise(
Richness = round(mean(Richness)),
.groups = "drop"
),
table(Species, Richness))
Richness
Species 0 1 2 3 4 5 6 7 8 9 10
1 80793 0 0 0 0 0 0 0 0 0 0
2 73294 1308 345 5846 0 0 0 0 0 0 0
3 60300 12419 913 7161 0 0 0 0 0 0 0
4 71970 7391 681 751 0 0 0 0 0 0 0
5 77639 2117 30 1007 0 0 0 0 0 0 0
6 80793 0 0 0 0 0 0 0 0 0 0
7 80793 0 0 0 0 0 0 0 0 0 0
8 80793 0 0 0 0 0 0 0 0 0 0
9 80793 0 0 0 0 0 0 0 0 0 0
10 62872 8776 1697 7448 0 0 0 0 0 0 0
11 80505 246 42 0 0 0 0 0 0 0 0
12 680 57 1 67 50 47 85 98 149 131 79428
13 80793 0 0 0 0 0 0 0 0 0 0
14 84 17 540 961 2714 3766 9410 3966 6741 12115 40479
15 71934 2740 85 6034 0 0 0 0 0 0 0
16 80793 0 0 0 0 0 0 0 0 0 0
17 47555 7152 340 5033 5534 8193 1609 1541 3836 0 0
18 80499 294 0 0 0 0 0 0 0 0 0
19 80686 107 0 0 0 0 0 0 0 0 0
20 80793 0 0 0 0 0 0 0 0 0 0
21 79313 1480 0 0 0 0 0 0 0 0 0
22 62 14 0 32 0 36 0 31 0 2888 77730
23 80793 0 0 0 0 0 0 0 0 0 0
24 80793 0 0 0 0 0 0 0 0 0 0
25 80793 0 0 0 0 0 0 0 0 0 0
26 80407 386 0 0 0 0 0 0 0 0 0
27 80793 0 0 0 0 0 0 0 0 0 0
28 80793 0 0 0 0 0 0 0 0 0 0
29 79711 114 387 581 0 0 0 0 0 0 0
30 23183 15781 13353 12433 16043 0 0 0 0 0 0
31 80793 0 0 0 0 0 0 0 0 0 0
32 4784 13308 14738 6522 8650 13321 6596 6339 2225 3634 676
33 34246 12708 12031 8137 10278 2444 722 227 0 0 0
34 80685 108 0 0 0 0 0 0 0 0 0
35 80600 193 0 0 0 0 0 0 0 0 0
36 56841 12740 11212 0 0 0 0 0 0 0 0
37 7020 40060 26246 7467 0 0 0 0 0 0 0
38 73593 7200 0 0 0 0 0 0 0 0 0
39 80793 0 0 0 0 0 0 0 0 0 0
40 78481 2312 0 0 0 0 0 0 0 0 0
41 80271 522 0 0 0 0 0 0 0 0 0
42 73553 7240 0 0 0 0 0 0 0 0 0
43 18796 19884 26680 11558 3875 0 0 0 0 0 0
44 30628 14552 5154 12331 12054 6074 0 0 0 0 0
45 77510 3283 0 0 0 0 0 0 0 0 0
46 70547 7937 160 2149 0 0 0 0 0 0 0
47 979 1092 10511 17784 13851 20944 11456 4176 0 0 0
48 52054 21980 6759 0 0 0 0 0 0 0 0
49 70981 9812 0 0 0 0 0 0 0 0 0
50 78655 2138 0 0 0 0 0 0 0 0 0
51 79858 935 0 0 0 0 0 0 0 0 0
52 69649 2937 894 5160 2153 0 0 0 0 0 0
53 79424 1369 0 0 0 0 0 0 0 0 0
54 78212 1554 1027 0 0 0 0 0 0 0 0
55 65177 15616 0 0 0 0 0 0 0 0 0
56 69587 11206 0 0 0 0 0 0 0 0 0
57 52153 20987 3657 3988 8 0 0 0 0 0 0
58 71255 9538 0 0 0 0 0 0 0 0 0
59 56304 24489 0 0 0 0 0 0 0 0 0
60 74416 6377 0 0 0 0 0 0 0 0 0
61 67723 12110 960 0 0 0 0 0 0 0 0
62 66212 4509 10072 0 0 0 0 0 0 0 0
63 75952 4841 0 0 0 0 0 0 0 0 0
64 76302 4491 0 0 0 0 0 0 0 0 0
65 80600 193 0 0 0 0 0 0 0 0 0
66 76287 4506 0 0 0 0 0 0 0 0 0
67 51275 22463 7055 0 0 0 0 0 0 0 0
68 26766 14846 7523 22838 7263 1557 0 0 0 0 0
69 70929 8209 1600 55 0 0 0 0 0 0 0
70 75050 5039 704 0 0 0 0 0 0 0 0
71 7364 4049 25788 11967 18316 13309 0 0 0 0 0
72 293 6015 31082 15185 8897 9826 6663 2832 0 0 0
73 65220 15573 0 0 0 0 0 0 0 0 0
74 76900 3893 0 0 0 0 0 0 0 0 0
75 78673 2120 0 0 0 0 0 0 0 0 0
76 60714 19313 662 104 0 0 0 0 0 0 0
77 66756 3870 5716 138 704 173 375 351 1496 1214 0
78 68931 11862 0 0 0 0 0 0 0 0 0
79 74538 6255 0 0 0 0 0 0 0 0 0
80 70951 9200 642 0 0 0 0 0 0 0 0
81 79192 1601 0 0 0 0 0 0 0 0 0
82 80444 349 0 0 0 0 0 0 0 0 0
83 79662 1131 0 0 0 0 0 0 0 0 0
84 73288 7505 0 0 0 0 0 0 0 0 0
85 78008 2785 0 0 0 0 0 0 0 0 0
86 8225 25572 17279 24066 5651 0 0 0 0 0 0
87 79538 1255 0 0 0 0 0 0 0 0 0
88 80538 255 0 0 0 0 0 0 0 0 0
89 38907 34291 7595 0 0 0 0 0 0 0 0
90 80793 0 0 0 0 0 0 0 0 0 0
[ reached getOption("max.print") -- omitted 10 rows ]
AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.06245 0.05166 0.04374 0.03644 0.03204 0.02918
sum(summary(PCA)$importance[2, ])
[1] 1.00007
ggplot2::autoplot(PCA, loadings = TRUE,
data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::mutate(
dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
)
PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)),
center = TRUE, scale. = TRUE, rank. = 25)
head(summary(PCAPA)$importance[2, ])
PC1 PC2 PC3 PC4 PC5 PC6
0.09404 0.07691 0.05775 0.05016 0.04514 0.03782
ggplot2::autoplot(PCAPA, loadings = TRUE,
data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
colour = "time")
sum(summary(PCA)$importance[2, ])
[1] 1.00007
ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)
Environments <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
retval <- data.frame(Environment = toString(i),
Time = time
)
retval <- cbind(retval, env)
colnames(retval) <- c("Environment", "Time",
paste0("Basal", 1:34),
paste0("Consumer", 35:100)
)
return(retval)
},
abund = result$Abundance %>% data.frame(
) %>% dplyr::mutate(
time = floor(time/100)*100
) %>% dplyr::group_by(
time
) %>% dplyr::summarise(
dplyr::across(.fns = ~ mean(.x))
),
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)
EnvironmentDistance <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard")
EnvDistClust <- hclust(EnvironmentDistance)
plot(EnvDistClust, labels = FALSE)
EnvironmentDistancePA <- Environments %>% dplyr::select(
-Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)
EnvDistClustPA <- hclust(EnvironmentDistancePA)
plot(EnvDistClustPA, labels = FALSE)